      subroutine alq4g2(coorr,coefr,
     & prmt,estif,emass,edamp,eload,num)
c .... coorr ---- nodal coordinate value
c .... coefr ---- nodal coef value
      implicit real*8 (a-h,o-z)
      dimension estif(12,12),elump(12),emass(12),
     & eload(12)
      dimension prmt(*),coef(3),coefr(4,3),coorr(2,4),coor(2)
      common /ralq4g2/ru(4,12),rv(4,12),rw(4,12),
     & cu(4,3),cv(4,3),cw(4,3)
c .... store shape functions and their partial derivatives
c .... for all integral points
      common /valq4g2/rctr(2,2),crtr(2,2),coefd(3,5),coefc(3,5)
      common /dalq4g2/ refc(2,4),gaus(4),
     & nnode,ngaus,ndisp,nrefc,ncoor,nvar,
     & nvard(3),kdord(3),kvord(12,3)
c .... nnode ---- the number of nodes
c .... nrefc ---- the number of numerical integral points
c .... ndisp ---- the number of unknown functions
c .... nrefc ---- the number of reference coordinates
c .... nvar ---- the number of unknown varibles var
c .... refc ---- reference coordinates at integral points
c .... gaus ---- weight number at integral points
c .... nvard ---- the number of var for each unknown
c .... kdord ---- the highest differential order for each unknown
c .... kvord ---- var number at integral points for each unknown
      fu=prmt(1)
      fv=prmt(2)
      fw=prmt(3)
      time=prmt(4)
      dt=prmt(5)
      imate=prmt(6)+0.5
      ielem=prmt(7)+0.5
      nelem=prmt(8)+0.5
      it=prmt(9)+0.5
      nmate=prmt(10)+0.5
      itime=prmt(11)+0.5
      ityp=prmt(12)+0.5
      if (num.eq.1) call alq4g2i
c .... initialize the basic data
      do 10 i=1,nvar
      eload(i)=0.0
      do 10 j=1,nvar
      estif(i,j)=0.0
10    continue
      do 999 igaus=1,ngaus
      call alq4g2t(nnode,nrefc,ncoor,refc(1,igaus),coor,coorr,
     & rctr,crtr,det,coefr)
c .... coordinate transfer from reference to original system
c .... rctr ---- Jacobi's matrix
c .... crtr ---- inverse matrix of Jacobi's matrix
      x=coor(1)
      y=coor(2)
      rx=refc(1,igaus)
      ry=refc(2,igaus)
      call ealq4g2(refc(1,igaus),coef,coorr,coefr,coefd)
c .... compute coef functions and their partial derivatives
      iu=(igaus-1)*3+1
      iv=(igaus-1)*3+1
      iw=(igaus-1)*3+1
      if (num.gt.1) goto 2
c .... the following is the shape function caculation
      call alq4g21(refc(1,igaus),ru(1,iu),rctr,crtr)
      call alq4g22(refc(1,igaus),rv(1,iv),rctr,crtr)
      call alq4g23(refc(1,igaus),rw(1,iw),rctr,crtr)
2     continue
c .... the following is the shape function transformation
c .... from reference coordinates to original coordinates
      call shapn(nrefc,ncoor,4,ru(1,iu),cu,crtr,1,3,3)
      call shapn(nrefc,ncoor,4,rv(1,iv),cv,crtr,1,3,3)
      call shapn(nrefc,ncoor,4,rw(1,iw),cw,crtr,1,3,3)
c .... the coef function transformation
c .... from reference coordinates to original coordinates
      call shapc(nrefc,ncoor,3,coefd,coefc,crtr,2,5,5)
      un=coef(1)
      vn=coef(2)
      wn=coef(3)
      weigh=det*gaus(igaus)
c .... the following is the stiffness computation
      do 202 i=1,4
      iv=kvord(i,1)
      do 201 j=1,4
      jv=kvord(j,1)
      stif=+cu(i,1)*cu(j,1)*0.0
      estif(iv,jv)=estif(iv,jv)+stif*weigh
201    continue
202    continue
c .... the following is the load vector computation
      do 501 i=1,4
      iv=kvord(i,1)
      stif=+cu(i,1)*fu
      eload(iv)=eload(iv)+stif*weigh
501   continue
      do 502 i=1,4
      iv=kvord(i,2)
      stif=+cv(i,1)*fv
      eload(iv)=eload(iv)+stif*weigh
502   continue
      do 503 i=1,4
      iv=kvord(i,3)
      stif=+cw(i,1)*fw
      eload(iv)=eload(iv)+stif*weigh
503   continue
999   continue
998   continue
      return
      end

      subroutine alq4g2i
      implicit real*8 (a-h,o-z)
      common /dalq4g2/ refc(2,4),gaus(4),
     & nnode,ngaus,ndisp,nrefc,ncoor,nvar,
     & nvard(3),kdord(3),kvord(12,3)
c .... initial data
c .... refc ---- reference coordinates at integral points
c .... gaus ---- weight number at integral points
c .... nvard ---- the number of var for each unknown
c .... kdord ---- the highest differential order for each unknown
c .... kvord ---- var number at integral points for each unknown
      ngaus=  4
      ndisp=  3
      nrefc=  2
      ncoor=  2
      nvar = 12
      nnode=  4
      kdord(1)=1
      nvard(1)=4
      kvord(1,1)=1
      kvord(2,1)=4
      kvord(3,1)=7
      kvord(4,1)=10
      kdord(2)=1
      nvard(2)=4
      kvord(1,2)=2
      kvord(2,2)=5
      kvord(3,2)=8
      kvord(4,2)=11
      kdord(3)=1
      nvard(3)=4
      kvord(1,3)=3
      kvord(2,3)=6
      kvord(3,3)=9
      kvord(4,3)=12
      refc(1,1)=5.773502692e-001
      refc(2,1)=5.773502692e-001
      gaus(1)=1.000000000e+000
      refc(1,2)=5.773502692e-001
      refc(2,2)=-5.773502692e-001
      gaus(2)=1.000000000e+000
      refc(1,3)=-5.773502692e-001
      refc(2,3)=5.773502692e-001
      gaus(3)=1.000000000e+000
      refc(1,4)=-5.773502692e-001
      refc(2,4)=-5.773502692e-001
      gaus(4)=1.000000000e+000
      end


      subroutine alq4g2t(nnode,nrefc,ncoor,refc,coor,coorr,
     & rc,cr,det,coefr)
      implicit real*8 (a-h,o-z)
      dimension refc(nrefc),rc(ncoor,nrefc),cr(nrefc,ncoor),a(5,10),
     & coorr(ncoor,nnode),coor(ncoor),coefr(nnode,*)
      call talq4g2(refc,coor,coorr,coefr,rc)
      n=nrefc
      m=n*2
      det = 1.0
      do 10 i=1,n
      do 10 j=1,n
      if (i.le.ncoor) a(i,j) = rc(i,j)
      if (i.gt.ncoor) a(i,j)=1.0
      a(i,n+j)=0.0
      if (i.eq.j) a(i,n+i) = 1.0
10    continue
c     write(*,*) 'a ='
c     do 21 i=1,n
c21   write(*,8) (a(i,j),j=1,m)
      do 400 i=1,n
      amax = 0.0
      l = 0
      do 50 j=i,n
      c = a(j,i)
      if (c.lt.0.0) c = -c
      if (c.le.amax) goto 50
      amax = c
      l = j
50    continue
      do 60 k=1,m
      c = a(l,k)
      a(l,k) = a(i,k)
      a(i,k) = c
60    continue
      c = a(i,i)
      det = c*det
      do 100 k=i+1,m
100   a(i,k) = a(i,k)/c
      do 300 j=1,n
      if (i.eq.j) goto 300
      do 200 k=i+1,m
200   a(j,k) = a(j,k)-a(i,k)*a(j,i)
c     write(*,*) 'i =',i,'  j =',j,'  a ='
c     do 11 ii=1,n
c11   write(*,8) (a(ii,jj),jj=1,m)
300   continue
400   continue
      do 500 i=1,nrefc
      do 500 j=1,ncoor
500   cr(i,j) = a(i,n+j)
c     write(*,*) 'a ='
c     do 22 i=1,n
c22   write(*,8) (a(i,j),j=1,m)
c     write(*,*) 'rc ='
c     do 24 i=1,ncoor
c24   write(*,8) (rc(i,j),j=1,nrefc)
c     write(*,*) 'cr ='
c     do 23 i=1,nrefc
c23   write(*,8) (cr(i,j),j=1,ncoor)
c     write(*,*) 'det =',det
      if (det.lt.0.0) det=-det
c     write(*,*) 'det =',det
8     format(1x,6f12.3)
      end

      subroutine alq4g21(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(2),shpr(4,3),rctr(2,2),crtr(2,2)
      external falq4g21
      rx=refc(1)
      ry=refc(2)
      call dshap(falq4g21,refc,shpr,2,4,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function falq4g21(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccalq4g2/ xa(4),ya(4),una(4),vna(4),
     &wna(4)
      common /valq4g2/ rctr(2,2),crtr(2,2),coefd(3,5),coefc(3,5)
      dimension refc(2)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      rx=refc(1)
      ry=refc(2)
      goto (1,2,3,4) n
1     falq4g21=+(+1.-rx)/2.*(+1.-ry)/2. 
      goto 1000
2     falq4g21=+(+1.+rx)/2.*(+1.-ry)/2. 
      goto 1000
3     falq4g21=+(+1.+rx)/2.*(+1.+ry)/2. 
      goto 1000
4     falq4g21=+(+1.-rx)/2.*(+1.+ry)/2. 
      goto 1000
1000  return
      end

      subroutine alq4g22(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(2),shpr(4,3),rctr(2,2),crtr(2,2)
      external falq4g22
      rx=refc(1)
      ry=refc(2)
      call dshap(falq4g22,refc,shpr,2,4,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function falq4g22(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccalq4g2/ xa(4),ya(4),una(4),vna(4),
     &wna(4)
      common /valq4g2/ rctr(2,2),crtr(2,2),coefd(3,5),coefc(3,5)
      dimension refc(2)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      rx=refc(1)
      ry=refc(2)
      goto (1,2,3,4) n
1     falq4g22=+(+1.-rx)/2.*(+1.-ry)/2. 
      goto 1000
2     falq4g22=+(+1.+rx)/2.*(+1.-ry)/2. 
      goto 1000
3     falq4g22=+(+1.+rx)/2.*(+1.+ry)/2. 
      goto 1000
4     falq4g22=+(+1.-rx)/2.*(+1.+ry)/2. 
      goto 1000
1000  return
      end

      subroutine alq4g23(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(2),shpr(4,3),rctr(2,2),crtr(2,2)
      external falq4g23
      rx=refc(1)
      ry=refc(2)
      call dshap(falq4g23,refc,shpr,2,4,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function falq4g23(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccalq4g2/ xa(4),ya(4),una(4),vna(4),
     &wna(4)
      common /valq4g2/ rctr(2,2),crtr(2,2),coefd(3,5),coefc(3,5)
      dimension refc(2)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      rx=refc(1)
      ry=refc(2)
      goto (1,2,3,4) n
1     falq4g23=+(+1.-rx)/2.*(+1.-ry)/2. 
      goto 1000
2     falq4g23=+(+1.+rx)/2.*(+1.-ry)/2. 
      goto 1000
3     falq4g23=+(+1.+rx)/2.*(+1.+ry)/2. 
      goto 1000
4     falq4g23=+(+1.-rx)/2.*(+1.+ry)/2. 
      goto 1000
1000  return
      end

      subroutine talq4g2(refc,coor,coorr,coefr,rc)
c .... compute coordinate value and Jacobi's matrix rc
c .... by reference coordinate value
      implicit real*8 (a-h,o-z)
      dimension refc(2),coor(2),coorr(2,4),coefr(4,3),rc(2,2)
      common /ccalq4g2/ x(4),y(4),un(4),vn(4),wn(4)
      external ftalq4g2
      do 100 n=1,4
      x(n)=coorr(1,n)
      y(n)=coorr(2,n)
100   continue
      do 200 n=1,4
      un(n)=coefr(n,1)
      vn(n)=coefr(n,2)
      wn(n)=coefr(n,3)
200   continue
      rx=refc(1)
      ry=refc(2)
      call dcoor(ftalq4g2,refc,coor,rc,2,2,1)
c .... coordinate value and their partial derivatives caculation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function ftalq4g2(refc,n)
c .... coordinate transfer function caculation
      implicit real*8 (a-h,o-z)
      dimension refc(2)
      common /ccalq4g2/ x(4),y(4),un(4),vn(4),wn(4)
      common /valq4g2/ rctr(2,2),crtr(2,2),coefd(3,5),coefc(3,5)
      rx=refc(1)
      ry=refc(2)
      goto (1,2) n
1     ftalq4g2=+(+(+1.-rx)/2.*(+1.-ry)/2.)*x(1)+(+(+1.+rx)/
     & 2.*(+1.-ry)/2.)*x(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*x(3)+(+
     & (+1.-rx)/2.*(+1.+ry)/2.)*x(4)
      goto 1000
2     ftalq4g2=+(+(+1.-rx)/2.*(+1.-ry)/2.)*y(1)+(+(+1.+rx)/
     & 2.*(+1.-ry)/2.)*y(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*y(3)+(+
     & (+1.-rx)/2.*(+1.+ry)/2.)*y(4)
      goto 1000
1000  return
      end

      subroutine ealq4g2(refc,coef,coorr,coefr,coefd)
c .... compute coef value and their partial derivatives
c .... by reference coordinate value
      implicit real*8 (a-h,o-z)
      dimension refc(2),coef(3),coorr(2,4),coefr(4,3),coefd(3,2)
      external fealq4g2
      rx=refc(1)
      ry=refc(2)
      call dcoef(fealq4g2,refc,coef,coefd,2,3,2)
c .... coef value and their partial derivatives caculation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function fealq4g2(refc,n)
c .... coef function caculation
      implicit real*8 (a-h,o-z)
      dimension refc(2)
      common /ccalq4g2/ xa(4),ya(4),un(4),vn(4),wn(4)
      common /valq4g2/ rctr(2,2),crtr(2,2),coefd(3,5),coefc(3,5)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      rx=refc(1)
      ry=refc(2)
      goto (1,2,3) n
1     fealq4g2=+(+(+1.-rx)/2.*(+1.-ry)/2.)*un(1)+(+(+1.+rx)/
     & 2.*(+1.-ry)/2.)*un(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*un(3)
     & +(+(+1.-rx)/2.*(+1.+ry)/2.)*un(4)
      goto 1000
2     fealq4g2=+(+(+1.-rx)/2.*(+1.-ry)/2.)*vn(1)+(+(+1.+rx)/
     & 2.*(+1.-ry)/2.)*vn(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*vn(3)
     & +(+(+1.-rx)/2.*(+1.+ry)/2.)*vn(4)
      goto 1000
3     fealq4g2=+(+(+1.-rx)/2.*(+1.-ry)/2.)*wn(1)+(+(+1.+rx)/
     & 2.*(+1.-ry)/2.)*wn(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*wn(3)
     & +(+(+1.-rx)/2.*(+1.+ry)/2.)*wn(4)
      goto 1000
1000  return
      end

